Tracking method for a radar system

ABSTRACT

A tracking method for a signal echo system, including generating a plurality of gates for respective propagation modes on the basis of a target state prediction for a dwell time, and generating a target state estimate for the dwell time on the basis of target measurement points which fall within the gates.

FIELD OF THE INVENTION

The present invention relates to a tracking method for a radar system,such as a phased array radar system or a bistatic radar system. Althoughthe following discusses use for radar systems, the invention could alsobe applied to other signal echo systems, such as sonar systems.

BACKGROUND OF THE INVENTION

Radar signals returned from a target allow information to be determinedconcerning the slant range, azimuth and speed of a target relative tothe receiving system of the radar system. The receiving system howevernormally receives a number of signals returned from the target whichhave different propagation paths or modes. Noise received by and inducedin the receiving system can also be mistaken for a return signal fromthe target and needs to be taken into account. Tracking methods havebeen employed which track a target on the basis of signals relating toone propagation mode. Yet selecting one propagation mode neglectsinformation relating to other modes which can be used to enhance theaccuracy and sensitivity of the tracking method.

BRIEF SUMMARY OF THE INVENTION

In accordance with the present invention there is provided a trackingmethod for a signal echo system, including:

generating a plurality of gates for respective propagation modes on thebasis of a target state prediction for a dwell time; and

generating a target state estimate for said dwell time on the basis oftarget measurement points which fall within said gates.

The present invention provides a tracking method for a signal echosystem, including:

obtaining target measurement points for a dwell time;

initiating tracking by obtaining an initial target state estimate fromat least one of said points;

determining a target state prediction for a subsequent dwell time on thebasis of said target state estimate;

generating a plurality of gates for respective propagation modes on thebasis of the target state prediction; and

generating a target state estimate for said subsequent dwell time on thebasis of target measurement points for said subsequent dwell time whichare within said gates.

The target state estimate may be generated by applying associationhypotheses to said measurement points in said gates and associationprobabilities to said hypotheses, obtaining conditional state estimatesfrom the measurement points for each hypothesis and summing saidconditional state estimates multiplied by said probabilities.

The tracking initiating step can advantageously be performed for aplurality of propagation modes to initiate a plurality of trackingfilters by generating a plurality of said target state estimates forsaid subsequent dwell time.

The present invention further provides a tracking method for a signalecho system, including extending a target state vector to includeadditional parameters associated with a plurality of propagation modes,and accounting for measurement uncertainty associated with propagationpath characteristics for said modes when updating target stateestimates.

BRIEF DESCRIPTION OF THE DRAWINGS

A preferred embodiment of the present invention is hereinafterdescribed, by way of example only, with reference to the accompanyingdrawings, wherein:

FIG. 1 is a schematic diagram of an over the horizon radar (OTHR)system, according to the prior art;

FIG. 2 is a diagram of a measurement frame of reference;

FIG. 3 is a diagram of possible propagation modes;

FIG. 4 is a graph of target tracks;

FIG. 5 is a diagram of measurement geometry of the system; and

FIG. 6 is a diagram of multiple gates for target state estimation.

DETAILED DESCRIPTION OF THE INVENTION

Bistatic radar systems employ separate transmitter and receiver sites,and include Over The Horizon Radar (OTHR) systems which directtransmission signals above the horizon for refraction by the ionosphere,known as skywave systems. OTHR systems also include surface wave radarsystems which propagate radar waves along the surface of saltwater, andrely on the receiving system being able to detect objects by the radarsignals reflected therefrom.

An OTHR system 2, as shown in FIG. 1, includes a receiving system 4 anda transmitting system 6. The transmitting system 6 comprises an in-linearray 7 of transmitting antennas located at the transmitter site and acontrol system 10 for feeding electrical signals to the antennas. Thereceiving system 4 comprises an in-line array 12 of receiving antennasand a control system 16 for processing the signals received by theantennas, which are located at the receiver site. OTHR systems includethe Jindalee Facility in Alice Springs (JFAS) and the U.S. Navy's ROTHRsystem.

The broad transmitting beam of the radar is directed towards areas ofthe ionosphere from which refracted signals are redirected to monitor atarget 3. The beam is effectively directed to a region or area in whicha target is located. A number of targets may be located in one regionand the receiver control system 16 is able to divide the energy returnedfrom the illuminated region into a dozen smaller beams which can theneach be divided into a plurality of range cells that are characterisedby a respective distance from the receiving system 4. This allows thereceiving system 4 to track a number of targets which are located in theilluminated region. The receive beams can also be divided into aplurality of velocity cells characterised by an object's velocityrelative to the receiving system 4. This allows targets to be separatedon the basis of their velocity if they cannot be separated on the basisof their distance from the receiving system 4. The transmitting andreceiving beams can be moved or swept in synchronism, through a numberof beam steer positions, with the time being spent at any given positionbeing referred to as the dwell time. Measurements obtained from theradar signals or echoes received during each dwell time are referred toas dwells.

The control software of the control system 16 is able to obtain fourparameters pertaining to a target from each dwell, and these are thepropagation path length or slant range (R), azimuth (A), Dopplerfrequency shift or radial speed (D) and signal strength based on asignal to noise ratio (SNR) measurement. These are referred to as theRAD or radar coordinates. The set of measurements from a dwell alsoincludes clutter and detections from other targets.

The dwells can be graphically represented by plotting them as candidatedetection points on a three dimensional axis, as shown in FIG. 2, fordwell t=k, where one axis represents R, the other A and the third the Dvalues. For any dwell time t=k of the order of 100 or 1000 candidatedetection points 50 may be determined by the receiving system 4. Some ofthe points 50 may correspond to a target and others may simply relate toclutter echoes or noise intrinsic in the transmitting system 6 orreceiving system 4. Clutter echoes arise from backscatter from theground or objects which are not of interest, such as meteors. The OTHRsystem 2 is also subject to multipath propagation in that there is morethan one single path for echoes returned from a target due to a numberof different ionospheric layers 54 at different heights 53 which refractechoes down to the receiving system 4. as shown in FIG. 3. There may beup to four different reflecting layers F₀, F₁, F₂ and F₃ resulting inseveral echoes returned from a target, corresponding to reflections fromcombinations of these layers. Propagation modes are described by thelayers from which the signal is refracted. For example, F₀-F₁ is thepropagation mode for a transmit path via layer F₀ and a receive path vialayer F₁, where T represents the target 3, as shown in FIG. 3. Whilstthe propagation path for a candidate detection point 50 is not known,the height of the different layers can be determined using commercialionospheric sounders which provides some information concerning therelationship between points of different propagation modes for the sametarget. Knowing the heights and properties of each layer gives anindication as to expected RAD measurements of different propagationmodes.

The state of the target, at a given dwell k, can be represented by$\begin{matrix}{{x(k)} = \begin{pmatrix}{r(k)} \\{\overset{.}{r}(k)} \\{a(k)} \\{\overset{.}{a}(k)}\end{pmatrix}} & (1)\end{matrix}$

where r is the range, a the azimuth, {dot over (r)} the range rate and{dot over (a)} the azimuth rate. Equations of motion can be used todescribe the target dynamics, for example, a constant velocity targetwould, if the time T between dwells were constant, obey

r(k)=r(0)+{dot over (r)}kT

a(k)=a(0)+{dot over (a)}kT.  (2)

This can be expressed in known state-space form as

x(k+1)=F(k)×(k)+v(k)  (3)

where F(k) is a known matrix, for instance in the case of a constantvelocity target $\begin{matrix}{{F(k)} = \begin{pmatrix}1 & T_{k} & 0 & 0 \\0 & 1 & 0 & 0 \\0 & 0 & 1 & T_{k} \\0 & 0 & 0 & 1\end{pmatrix}} & (4)\end{matrix}$

where T_(k) is the time between dwells k and k+1. The term ν(k)represents zero-mean, white Gaussian process noise as used in standardKalman filtering. The covariance matrix Q(k) of ν(k) is assumed to beknown.

A currently used tracking method, based on the probabilistic dataassociation (PDA) filter, as described in Y. Bar-Shalom and T. E.Fortmann, “Tracking and Data Association”, Academic Press, 1988.performs tracking in the radar coordinates R, D, A, A, as illustrated inFIG. 2. A track is initiated by selection of a single, noisy measurement50 with the unknown azimuth rate A being initially set to correspondwith a hypothesised azimuth crossing rate, usually zero. Furthermeasurement selection is accomplished by taking only those measurementswhich fall inside a validation gate 70 around the next expected positionof the target measurement. This method does not require knowledge of themapping of radar to ground coordinates during tracking. A disadvantageof this method is that it fails to use the information conveyed bymultiple detections arising from multipath propagation. Also thepresence of multipath propagation may cause multiple tracks 60, 62 and64 to be generated for a single target, as shown in FIG. 4, whentracking is performed in radar coordinates using conventional filterssuch as PDA. If the tracks 62 and 64 closely conform with the expectedseparations for the hypothesised modes, they can be considered to relateto the same target 52 of FIG. 5, whereas a track 66 which divergesexcessively can be dismissed as corresponding to another target or toclutter. Such a situation, commonly arising in conventional PDAtracking, requires a fusion or clustering operation to group multimodetracks pertaining to the same target together. This allows a track to beidentified with a particular propagation mode. A further stage ofcoordinate registration is then required to map the tracks to groundcoordinates for geographical display to the radar operators.

The preferred embodiment described herein uses explicit knowledge of theionospheric structure including virtual heights, as provided byionospheric sounders or by other means, to account for and takeadvantage of multipath propagation during track initiation and tracking.This is distinct from conventional approaches which only expect a singledetection per target and are unable to benefit from the additionaltarget-related information conveyed by multipath detections. The gain intracking performance arising from multipath detections of a singletarget is important when the probability of target detection via some orall of the various propagation modes is low.

The target state is taken to be as in equation (1), where r is theground range 8 across the surface of the earth 9, a is the true azimuth,{dot over (r)} is the ground range rate and a is the true azimuth rate.The true azimuth a is the complement of the angle {overscore (a)}, whichis the angle between the projected ground range r and the axis of thereceiver array 12, as shown in FIG. 5, i.e. a=(90°−{overscore (a)}).Tracking is performed in ground coordinates, although other frames ofreference, for instance a preferred propagation mode, may be used todescribe the target dynamics and relate these to the other measurementcoordinates.

The conversion between the ground and radar coordinates can berepresented as $\begin{matrix}{\begin{pmatrix}{R(k)} \\{A(k)} \\{D(k)}\end{pmatrix} = {H\left( {{r(k)},{a(k)},{{\overset{.}{r}(k)};h_{r}},\quad h_{t}} \right)}} & (5)\end{matrix}$

where, at time k, R is the measured slant range, A the measured azimuth,D the Doppler speed (slant range rate), h_(r) the virtual ionosphericheight 53 on the receive path, and h_(t) the virtual ionospheric height54 on the transmit path, as shown in FIG. 5. The slant range R may bedefined as one half of the total path length from the transmitter 7 viathe target 52 to the receiver 4. The measured azimuth or coning angle Ais the complement of the angle {overscore (A)} between the incoming ray57 and the axis of receiver array axis 12. The Doppler speed D isproportional to the rate of change of the total path length.

The various propagation modes can be labelled according to thecorresponding outbound and return propagation mode combination F₀-F₀,F₀-F₁, F₁-F₀, . . . , F₂-F₂ for a target 52, as shown in FIG. 3. Forfour possible ionospheric layers F₀, F₁, F₂, F₃ with heights h₀, h₁, h₂,h₃, these modes may be numbered from 1 to 16 respectively. Hence we maywrite the measurement process for the various propagation modes in termsof the target state x(k) as $\begin{matrix}{{y(k)} = \left\{ \begin{matrix}{{H_{1}\left( {x(k)} \right)} + {{w_{1}(k)}\quad \text{for mode}\quad F_{0}} - F_{0}} \\{{H_{2}\left( {x(k)} \right)} + {{w_{2}(k)}\quad \text{for mode}\quad F_{0}} - F_{1}} \\\vdots \\{{H_{n}\left( {x(k)} \right)} + {{w_{n}(k)}\quad \text{for mode}\quad n}}\end{matrix} \right.} & (6)\end{matrix}$

where H₁(x(k))=H(r(k), a(k), {dot over (r)}(k); h₀, h₀),H₂(x(k))=H(r(k), a(k), {dot over (r)}(k); h₀, h₁), etc., and the assumednumber of possible propagation modes n may vary with time. In the above,w_(i)(k) is a zero-mean, white Gaussian sequence with known covarianceR_(i)(k) representing the assumed measurement noise terms. The actualform of the non-linear measurement functions H_(i)(·) above isdetermined by the geometry of the ionospheric model as shown in FIG. 5,and will depend on the virtual heights of the ionospheric layers h_(r)and h_(t) 53 and 55, and the location and separation of the receiver andtransmitter arrays 7 and 12, among other factors.

Since the virtual ionospheric heights h_(i) in FIG. 5 may only beapproximately known, but are assumed to vary slowly in comparison withthe target dynamics, they can be included in the state vector x(k) andestimated along with the dynamical variables describing the target. Inthis case we have instead of equation (1) $\begin{matrix}{{x(k)} = \begin{pmatrix}{r(k)} \\{\overset{.}{r}(k)} \\{a(k)} \\{\overset{.}{a}(k)} \\{h_{1}(k)} \\\vdots \\{h_{n}(k)}\end{pmatrix}} & (7)\end{matrix}$

with each virtual height satisfying an equation of the form

h_(i)(k+1)=h_(i)(k)+ν_(i)(k)  (8)

where ν_(i)(k) is a small process noise term.

Converting to the ground frame of reference requires the selection of anoutbound and return propagation mode combination F_(t) and F_(r) withcorresponding virtual heights h_(t) and h_(r).

The inverse transformation to equation (5) can be represented by$\begin{matrix}{\begin{pmatrix}{r(k)} \\{a(k)} \\{\overset{.}{r}(k)}\end{pmatrix} = {H^{- 1}\left( {{R(k)},{A(k)},{{D(k)};h_{r}},h_{t}} \right)}} & (9)\end{matrix}$

and follows from the assumed geometry indicated in FIG. 5.

Hereinafter the state prediction and associated prediction covarianceare denoted by {circumflex over (x)}(k|k−1) and P(k|k−1) and an updatedstate estimate and state error covariance are denoted by {circumflexover (x)}(k|k) and P(k|k).

At some arbitrary time 0, tracking is initiated by selecting an initialpoint 50 which may correspond to a hitherto unobserved target. Since thepropagation mode which gave rise to this measurement is a prioriunknown, an initial target state estimate {circumflex over (x)}(0|0) forequation (1) cannot be inferred from equation (9) unless a givenpropagation mode, or equivalently the ionospheric heights for thetransmit and receive paths, is assumed. The preferred method istherefore to initialise n tracking filters, one for each possibleinitial propagation mode. Each filter assumes a particular initialpropagation mode with corresponding virtual heights h_(r) and h_(t) inorder to assign its initial state estimate using equation (9) based onthe first measurement point 50. The estimate of the initial targetazimuth rate is set to some starting value, usually zero. An initialstate error covariance P(0|0) is also assigned and is taken to be largeenough to cover the initial uncertainty in target position and velocity.Other methods of initialisation are possible using data from more than asingle radar dwell; but the previously described method is the simplestamong these. Of the n filters initiated from the measurement 50, thefilter based on the correct initial propagation mode assumption can beexpected to perform the best and thus its state estimates would be moreaccurate (in the sense of having smaller errors on average) than thoseof the other filters initiated with it. As the processing proceeds, itbecomes clear by observation of the state estimates which, if any, ofthe n filters initiated as above is compatible with a target whosedynamical model is assumed to be as expressed in equation (2).

The recursive processing required by each tracking filter, initialisedas above, is now described. The aim of the processing is to compute, ina recursive manner, approximate conditional mean {circumflex over(x)}(k|k) and covariance P(k|k) estimates of the target state, based onthe measurement data, including virtual ionospheric height measurements,up to time k, Y(1), . . . , Y(k), where Y(i) represents the set ofmeasurements received in dwell i. The estimated target track is providedby plotting the range and azimuth values from {circumflex over(x)}(k|k). The accuracy of the track is indicated by the size of thestandard deviations which can be obtained from the state errorcovariance P(k|k).

The dynamical target model of equation (3) is used to predict where eachmeasurement would appear during the next dwell in the absence ofmeasurement noise under each propagation mode. The state prediction{circumflex over (x)}(1|0) at time 1 is given, in the usual manner ofKalman filtering, as described in Y. Bar-Shalom and T. E. Fortmann,“Tracking and Data Association”, Academic Press, 1988, as

{circumflex over (x)}(1|0)=F(0){circumflex over (x)}(0|0)  (10)

with associated covariance

P(1|0)=F(0)P(0|0)F′(0)+Q(0)  (11)

where F′ is the transpose of the transition matrix F in equation (3).Instead of generating one gate 70, {circumflex over (x)}(1|0) is used togenerate n gates 72, 74 and 76 in the measurement space for eachtracking filter, corresponding to the respective propagation modesF₀-F₀, F₀-F₁, . . . , etc., as shown in FIG. 6. The measurementpredictions for the respective propagation modes are therefore$\begin{matrix}\begin{matrix}{{{\hat{y}}_{1}\left( 1 \middle| 0 \right)} = {H_{1}\left( {\hat{x}\left( 1 \middle| 0 \right)} \right)}} \\{{{\hat{y}}_{2}\left( 1 \middle| 0 \right)} = {H_{2}\left( {\hat{x}\left( 1 \middle| 0 \right)} \right)}} \\\vdots \\{{{\hat{y}}_{n}\left( 1 \middle| 0 \right)} = {{H_{n}\left( {\hat{x}\left( 1 \middle| 0 \right)} \right)}.}}\end{matrix} & (12)\end{matrix}$

The associated measurement prediction covariances are $\begin{matrix}\begin{matrix}{{S_{1}(1)} = {{{J_{1}(1)}{P\left( 1 \middle| 0 \right)}{J_{1}^{\prime}(1)}} + {R_{1}(1)}}} \\{{S_{2}(1)} = {{{J_{2}(1)}{P\left( 1 \middle| 0 \right)}{J_{2}^{\prime}(1)}} + {R_{2}(1)}}} \\\vdots \\{{S_{n}(1)} = {{{J_{n}(1)}{P\left( 1 \middle| 0 \right)}{J_{n}^{\prime}(1)}} + {R_{n}(1)}}}\end{matrix} & (13)\end{matrix}$

where J_(i)(1) is the Jacobian matrix of the non-linear measurementfunction H_(i)(·) in equation (6) evaluated at the state predictions{circumflex over (x)}(1|0). The validation gate for each propagationmode is an ellipsoidal region in RAD space defined by

G_(i)(1)={yεIR³:[y−ŷ_(i)(1|0)]′S_(i)(1)⁻¹[y−ŷ_(i)(1|0]<γ_(i)}  (14)

where γ_(i), defines the size of the validation gate. The probabilitythat a target falls inside the gate i is denoted P_(G) ^(i), while theprobability of detecting a target via the ith propagation mode isdenoted P_(D) ^(i). This is illustrated in FIG. 6 for three gates 72, 74and 76 centred on measurement predictions 82, 84 and 86 for threepropagation modes F₀-F₀, F₀-F₁, . . . , etc. The gates may or may notoverlap. The validation region is defined as the union of the validationgates or some region which includes their union. Points which fallinside the validation gates are accepted as possibly relating to thetarget 52 and are used together with the state prediction {circumflexover (x)}(1|0) in order to update the state estimate {circumflex over(x)}(0|0) to yield {circumflex over (x)}(1|1). The corresponding stateerror covariance is also updated to P(1|1). This process is recursiveand can be represented as follows:

{circumflex over (x)}(0|0)→{circumflex over (x)}(1|0)→{circumflex over(x)}(1|1)→{circumflex over (x)}(2|1)→ . . .

P(0|0)→P(1|0)→P(1|1)→P(2|1)→ . . .  (15)

The state estimate {circumflex over (x)}(k|k) is an approximateminimum-mean-square-error estimate of the target state x(k) based on allthe information 48 from dwells 0 through k of the form given in FIG. 2including multiple detections of the same target due to multipathpropagation. The estimate is approximate because it assumes that theprobability density function of the true target state is Gaussianconditioned on all the measurement data.

To determine the updated target state {circumflex over (x)}(k|k) and itscovariance P(k|k), the measurements falling within the gates 72, 74, 76.etc. are used in a probabilistic data association framework as describedin Y. Bar-Shalom and T. E. Fortmann, “Tracking and Data Association”,Academic Press, 1988. which, in addition to consideration of ameasurement being from a target or due to clutter, includes associationhypotheses for the possible propagation modes which may have producedthe measurements. A target existence or confidence model is alsoincorporated in the filter as described in S. B. Colegrove, A. W. Davisand J. K. Ayliffe, “Track Initiation and Nearest Neighbours Incorporatedinto Probabilistic Data Association”, J. Elec. and Electronic Eng.,Australia, Vol. 6. No. 3, pp. 191-198, 1986, to aid in track maintenance(confirmation, deletion, etc.). The probability that the target existsat time k given data to time k is denoted P_(E)(k|k). Target existenceis modelled as a 2-state Markov Chain so that the predicted probabilityof target existence P_(E)(k|k−1) satisfies

P_(E)(k|k−1)=Δ₀P_(E)(k−1|k−1)+Δ₁{1−P_(E)(k−1|k−1)}  (16)

where the two transition probabilities Δ₀ and Δ₁ are defined by

Δ₀=Pr(target exists at time k|target exists at time k−1)

Δ₁=Pr(target exists at time k|target does not exist at time k−1)

An arbitrary initial value of P_(E)(0|0)=0.5 is assumed.

As an illustration of the filtering procedure, consider gates 72 and 74associated with propagation modes F₀-F₀ and F₀-F₁, and gate 76associated with propagation mode F₁-F₁, with respective centres given bythe measurement predictions ŷ₁(k|k−1), ŷ₂(k|k−1) and y₃(k|k−1), 82, 84and 86. We will number these propagation modes as 1, 2 and 3,respectively when referring to the measurement predictions. Suppose thatthe gate 72 contains two measurements y₁, Y₂ 90 and that gate 74contains one measurement y₃ 92, while gate 76 does not contain anymeasurements. The 7 association hypothesis (numbered from −1 to 5) whichcan be applied are:

(−1) The target does not exist.

(0) The target exists but all validated measurements y₁, y₂ and y₃ areclutter.

(1) y₃ and y₃ are clutter, y₂ is a target detection via propagation modeF₀-F₀.

(2) y₂ and y₃ are clutter, y₁ is a target detection via propagation modeF₀-F₀.

(3) y₃ is a target detection via F₀-F₁, and both y₁, and y₂ are clutter.

(4) y₃ is a target detection via F₀-F₁, y₁ is a target detection viaF₀-F₀ and y₂ is clutter.

(5) y₃ is a target detection via F₀-F₁, Y₂ is a target detection viaF₀-F₀ and y₁ is clutter.

For each of the possible associated hypothesis above, a conditionaltarget state estimate {circumflex over (x)}_(i)(k|k) can be formed fromthe predicted state estimate {circumflex over (x)}(k|k−1) using theextended Kalman filter theory, as described in G. W. Pulford and R. J.Evans, “Probabilistic Data Association for Systems with MultipleSimultaneous Measurements”, Automatica, Vol. 32. No. 9. pp. 1311-1316,1996. Omitting some time indexes and writing {overscore (x)}={circumflexover (x)}(k|k−1) for equation (10), {overscore (P)}=P(k|k−1) forequation (11), and {overscore (y)}_(i)=ŷ_(i)(k|k−1) for equation (12),the conditional state estimates in this case are given by$\begin{matrix}{{{\hat{x}}_{- 1}\left( k \middle| k \right)} = \overset{\_}{x}} & (17) \\{{{\hat{x}}_{0}\left( k \middle| k \right)} = \overset{\_}{x}} & \quad \\{{{\hat{x}}_{1}\left( k \middle| k \right)} = {\overset{\_}{x} + {\overset{\_}{P}J_{1}^{\prime}S_{1}^{- 1}\left\{ {y_{2} - {\overset{\_}{y}}_{1}} \right\}}}} & \quad \\{{{\hat{x}}_{2}\left( k \middle| k \right)} = {\overset{\_}{x} + {\overset{\_}{P}J_{1}^{\prime}S_{1}^{- 1}\left\{ {y_{1} - {\overset{\_}{y}}_{1}} \right\}}}} & \quad \\{{{\hat{x}}_{3}\left( k \middle| k \right)} = {\overset{\_}{x} + {\overset{\_}{P}J_{2}^{\prime}S_{2}^{- 1}\left\{ {y_{3} - {\overset{\_}{y}}_{2}} \right\}}}} & \quad \\{{{\hat{x}}_{4}\left( k \middle| k \right)} = {\overset{\_}{x} + {{\overset{\_}{P}\left( {J_{1}^{\prime}J_{2}^{\prime}} \right)}\begin{pmatrix}S_{1} & {J_{1}\overset{\_}{P}J_{2}^{\prime}} \\{J_{2}\overset{\_}{P}J_{1}^{\prime}} & S_{2}\end{pmatrix}^{- 1}\begin{pmatrix}y_{1} & {- {\overset{\_}{y}}_{1}} \\y_{3} & {- {\overset{\_}{y}}_{2}}\end{pmatrix}}}} & \quad \\{{{\hat{x}}_{5}\left( k \middle| k \right)} = {\overset{\_}{x} + {{\overset{\_}{P}\left( {J_{1}^{\prime}J_{2}^{\prime}} \right)}\begin{pmatrix}S_{1} & {J_{1}\overset{\_}{P}J_{2}^{\prime}} \\{J_{2}\overset{\_}{P}J_{1}^{\prime}} & S_{2}\end{pmatrix}^{- 1}\begin{pmatrix}y_{2} & {- {\overset{\_}{y}}_{1}} \\y_{3} & {- {\overset{\_}{y}}_{2}}\end{pmatrix}}}} & \quad\end{matrix}$

where the terms J_(i); and S_(i) are as in equation (13). Thecorresponding conditional state error covariances P_(i)(k|k) are givenby $\begin{matrix}{{{P_{- 1}\left( k \middle| k \right)} = {c\overset{\_}{P}}}{{P_{0}\left( k \middle| k \right)} = \overset{\_}{P}}{{P_{1}\left( k \middle| k \right)} = {\overset{\_}{P} - {\overset{\_}{P}J_{1}^{\prime}S_{1}^{- 1}J_{1}\overset{\_}{P}}}}{{P_{2}\left( k \middle| k \right)} = {\overset{\_}{P} - {\overset{\_}{P}J_{1}^{\prime}S_{1}^{- 1}J_{1}\overset{\_}{P}}}}{{P_{3}\left( k \middle| k \right)} = {\overset{\_}{P} - {\overset{\_}{P}J_{2}^{\prime}S_{2}^{- 1}J_{2}\overset{\_}{P}}}}{{P_{4}\left( k \middle| k \right)} = {\overset{\_}{P} - {{\overset{\_}{P}\left( {J_{1}^{\prime}J_{2}^{\prime}} \right)}\begin{pmatrix}S_{1} & {J_{1}\overset{\_}{P}J_{2}^{\prime}} \\{J_{2}\overset{\_}{P}J_{1}^{\prime}} & S_{2}\end{pmatrix}\begin{pmatrix}J_{1} \\J_{2}\end{pmatrix}\overset{\_}{P}}}}{{P_{5}\left( k \middle| k \right)} = {\overset{\_}{P} - {{\overset{\_}{P}\left( {J_{1}^{\prime}J_{2}^{\prime}} \right)}\begin{pmatrix}S_{1} & {J_{1}\overset{\_}{P}J_{2}^{\prime}} \\{J_{2}\overset{\_}{P}J_{1}^{\prime}} & S_{2}\end{pmatrix}\begin{pmatrix}J_{1} \\J_{2}\end{pmatrix}\overset{\_}{P}}}}} & (18)\end{matrix}$

where c≧1 is a scaling factor reflecting the increased uncertainty inthe case that the target does not exist.

The computation of the probability of each associated hypothesis, calledthe association probability, can be illustrated by assuming uniformlydistributed clutter measurements in the radar measurement space, and aPoisson model, as described in Y. Bar-Shalom and T. E. Fortmann,“Tracking and Data Association”, Academic Press, 1988, with spatialdensity λ for the number of clutter points inside the validation region.Also, we let the probability of target detection P_(D) via anypropagation mode be identical, and the gate probabilities P_(G) beidentical. P_(D) and P_(G) are parameters that are given values whichare selected to extract optimum performance from the tracking methodgiven the operating characteristics of the system 2. The total volume ofthe validation gates at time k is V_(k). If two or more of the gatesoverlap, V_(k) can be approximated, for instance, as the volume of thelargest gate. The association probabilities β_(i)(k), defined as theprobability of the respective association hypothesis i conditioned onall measurement data up to the current time k, can be expressed asdescribed in G. W. Pulford and R. J. Evans, “Probabilistic DataAssociation for Systems with Multiple Simultaneous Measurements”,Automatica, Vol. 32. No. 9. pp. 1311-1316, 1996, by: $\begin{matrix}\begin{matrix}{{\beta_{- 1}(k)} = \quad {{\delta^{- 1}(k)}\left\{ {1 - {{P_{E}\left( {k\left. {k - 1} \right)} \right\}}{\exp \left( {{- \lambda}\quad V_{k}} \right)}{\lambda^{3}/{3!}}}} \right.}} \\{{\beta_{0}(k)} = \quad {{\delta^{- 1}(k)}\left( {1 - {P_{D}P_{G}}} \right)^{3}\lambda^{3}{\exp \left( {{- \lambda}\quad V_{k}} \right)}{{P_{E}\left( {k{{k - 1}}} \right)}/{3!}}}} \\{{\beta_{1}(k)} = \quad {{\delta^{- 1}(k)}{P_{D}\left( {1 - {P_{D}P_{G}}} \right)}^{2}\lambda^{2}{\exp \left( {{- \lambda}\quad V_{k}} \right)}}} \\{\quad {N\left\{ {{y_{2};{\overset{\_}{y}}_{1}},S_{1}} \right\} {{P_{E}\left( {k{{k - 1}}} \right)}/\left( {3 \times {2!}} \right)}}} \\{{\beta_{2}(k)} = \quad {{\delta^{- 1}(k)}{P_{D}\left( {1 - {P_{D}P_{G}}} \right)}^{2}\lambda^{2}{\exp \left( {{- \lambda}\quad V_{k}} \right)}}} \\{\quad {N\left\{ {{y_{1};{\overset{\_}{y}}_{1}},S_{1}} \right\} {{P_{E}\left( {k{{k - 1}}} \right)}/\left( {3 \times {2!}} \right)}}} \\{{\beta_{3}(k)} = \quad {{\delta^{- 1}(k)}{P_{D}\left( {1 - {P_{D}P_{G}}} \right)}^{2}\lambda^{2}{\exp \left( {{- \lambda}\quad V_{k}} \right)}}} \\{\quad {N\left\{ {{y_{3};{\overset{\_}{y}}_{2}},S_{2}} \right\} {{P_{E}\left( {k{{k - 1}}} \right)}/\left( {3 \times {2!}} \right)}}} \\{{\beta_{4}(k)} = \quad {{\delta^{- 1}(k)}{P_{D}^{2}\left( {1 - {P_{D}P_{G}}} \right)}{{\lambda exp}\left( {{- \lambda}\quad V_{k}} \right)}}} \\{\quad {N\left\{ {{y_{1};{\overset{\_}{y}}_{1}},S_{1}} \right\} N\left\{ {{y_{3};{\overset{\_}{y}}_{2}},S_{2}} \right\} {{P_{E}\left( {k{{k - 1}}} \right)}/\left( {2 \times {1!}} \right)}}} \\{{\beta_{5}(k)} = \quad {{\delta^{- 1}(k)}{P_{D}^{2}\left( {1 - {P_{D}P_{G}}} \right)}{{\lambda exp}\left( {{- \lambda}\quad V_{k}} \right)}}} \\{\quad {N\left\{ {{y_{2};{\overset{\_}{y}}_{1}},S_{1}} \right\} N\left\{ {{y_{3};{\overset{\_}{y}}_{2}},S_{2}} \right\} {{P_{E}\left( {k{{k - 1}}} \right)}/\left( {2 \times {1!}} \right)}}}\end{matrix} & (19)\end{matrix}$

where N{y; {overscore (y)}, S} is a multivariate Gaussian density inywith mean {overscore (y)} and covariance S, and δ(k) is a normalisationconstant, chosen to ensure that the association probabilities sum tounity.

The updated target state estimate for the filter is obtained by summingthe conditional state estimates with weightings determined by theirrespective association probabilities as $\begin{matrix}{\hat{x}\left( {{k\left. k \right)} = {\left\{ {{\beta_{- 1}(k)} + {\beta_{0}(k)}} \right\} {\hat{x}\left( {{k\left. {k - 1} \right)} + {\sum\limits_{i = 1}^{5}{{\beta_{i}(k)}{{\hat{x}}_{i}\left( {k{\left. k \right).}} \right.}}}} \right.}}} \right.} & (20)\end{matrix}$

The state error covariance P(k|k) is obtained using standard techniquesfrom Gaussian mixtures described in Y. Bar-Shalom and T. E. Fortmann,“Tracking and Data Association”, Academic Press, 1988 as $\begin{matrix}\begin{matrix}{{P\left( {k{k}} \right)} = \quad {\left\{ {{c\quad {\beta_{- 1}(k)}} + {\beta_{0}(k)}} \right\} {P\left( {{k\left. {k - 1} \right)} +} \right.}}} \\{\quad {\left\{ {{\beta_{- 1}(k)} + {\beta_{0}(k)}} \right\} {\hat{x}\left( {k\left. {k - 1} \right){\hat{x}\left( {{k\left. {k - 1} \right)^{\prime}} - {\hat{x}\left( {k\left. k \right){\hat{x}\left( {{k\left. k \right)^{\prime}} +} \right.}} \right.}} \right.}} \right.}}} \\{\quad {\sum\limits_{i = 1}^{5}{{\beta_{i}(k)}\left\{ {P_{i}\left( {{k\left. k \right)} + {{\hat{x}}_{i}\left( {k\left. k \right){{{\hat{x}}_{i}\left( {k\left. k \right)^{\prime}} \right\}}.}} \right.}} \right.} \right.}}}\end{matrix} & (21)\end{matrix}$

The updated probability of target existence P_(E)(k|k) is obtained as

P_(E)(k|k)=1−β⁻¹(k).  (22)

Track maintenance is achieved by thresholding the target existenceprobability according to $\begin{matrix}{P_{E}\left( {\left. {k{k}} \right) < P_{DEL}}\Rightarrow{{delete}\quad {track}P_{E}\left( {\left. {k{k}} \right) > P_{CON}}\Rightarrow{{{confirm}\quad {track}P_{DEL}} \leq {P_{E}\left( {\left. {k{k}} \right) \leq P_{CON}}\Rightarrow{{retain}\quad {track}\quad {as}\quad {tentative}} \right.}} \right.} \right.} & (23)\end{matrix}$

where P_(DEL) and P_(CON) are the track maintenance thresholds. NoteP_(DEL)<P_(CON). Since P_(E)(k|k) may vary considerably from dwell todwell, it is better to use the average value of P_(E)(k|k) over the lastfew dwells for track confirmation in equation (23).

The above method is easily extended to arbitrary numbers of measurementsfalling inside the validation gates and to arbitrary numbers ofpropagation modes. Arbitrary clutter probability density functions andnon-identical gate and detection probabilities are also easilyaccommodated within this framework.

A tracking filter as described above has been implemented in softwareusing the C programming language and executed on a Digital EquipmentCorporation 175 MHz Alpha workstation. The preferred implementationassumes 4 propagation modes corresponding to F₀-F₀, F₀-F₁, F₁-F₀ andF₁-F₁. The virtual heights of the F₀ and F₁. ionospheric layers areincluded as state variables in equation (7) and these are estimated fromnoisy measurements along with the range, azimuth, range rate and azimuthrate of the target.

Many modifications will be apparent to those skilled in the art withoutdeparting from the scope of the present invention as herein describedwith reference to the accompanying drawings.

What is claimed is:
 1. A tracking method for a signal echo system,including: generating a plurality of gates for respective propagationmodes on the basis of a target state prediction for a dwell time; andgenerating a target state estimate for said dwell time on the basis oftarget measurement points which fall within said gates.
 2. A trackingmethod as claimed in claim 1, including: obtaining initial targetmeasurement points for an initial dwell time; initiating tracking byobtaining an initial target state estimate from at least one of saidinitial points; and determining said target state prediction for a dwelltime subsequent to said initial time on the basis of said initial targetstate estimate.
 3. A tracking method as claimed in claim 2, wherein atarget state estimate is generated by applying association hypotheses tosaid measurement points in said gates and association probabilities tosaid hypotheses, obtaining conditional state estimates from themeasurement points for each hypothesis and summing said conditionalstate estimates multiplied by said probabilities.
 4. A tracking methodas claimed in claim 3, wherein the probability of existence P_(E) of atarget track is obtained from at least one of said associationprobabilities and if P_(E) is less than a predetermined threshold atarget track maintained using said target state estimate is deleted. 5.A tracking method as claimed in claim 4, wherein said measurement pointsare candidate detection points in RAD space obtained from dwells.
 6. Atracking method as claimed in claim 5, wherein said target stateprediction is obtained from said target state estimate using linearequations of motion.
 7. A tracking method as claimed in claim 6, whereinsaid gates are validation gates having an ellipsoidal shape in RAD spaceand obtained by transposing said target state prediction to RAD spacefor respective propagation modes to obtain measurement predictions andassociated prediction covariances for respective propagation modesdefining said validation gates.
 8. A tracking method as claimed in claim7, wherein said hypotheses include target does not exist, themeasurement points in said gates represent clutter, and a measurementpoint in at least one of said gates represents a target detection.
 9. Atracking method as claimed in claim 2, wherein the tracking initiatingstep is performed for a plurality of propagation modes to initiate aplurality of tracking filters by generating a plurality of said initialtarget state estimates.
 10. A tracking method for a signal echo system,including extending a target state vector to include additionalparameters associated with a plurality of propagation modes, generatinga plurality of gates for said propagation modes, and accounting formeasurement uncertainty associated with propagation path characteristicsfor said modes when updating target state estimates for a dwell time onthe basis of measurements which fall within the gates.